%% Set paths:

clear;clc;close all

%% Read data

data_orig           = readtable('bhc_and_commercial_hh_min10_branches_for_matlab_nobranchconstraints.csv');
data_orig.z         = data_orig.assets./data_orig.brnumber_hh; % Create additional variables

%% Figures OA.2.:Figures: CDF for each year post DF
amin = 3;
amax = 250;
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = data(data.post == 0,:).assets;
apost = data(data.post == 2,:).assets;


data_post = data_orig(data_orig.post == 2, :);
data_post.year = floor(data_post.datequarter_num);
data_post = data_post(data_post.year >=2011, :);

yearlist = 2011:2018;

apre_9_11_struc = {};
apost_9_11_struc = {};
apost_48_52_struc = {};
x_pre_struc = {};
x_post_struc = {};
x_post_struct_50b = {};


for y = 1:length(yearlist)
    apost_9_11_struc{y} = data_post(data_post.assets <=11 & data_post.assets>= 9 & data_post.year == yearlist(y),:).assets;
    apost_48_52_struc{y} = data_post(data_post.assets <=52 & data_post.assets>= 48 & data_post.year == yearlist(y),:).assets;

    x_post_struc{y}     = (1:length( apost_9_11_struc{y}))/length( apost_9_11_struc{y});
    x_post_struct_50b{y} = (1:length(  apost_48_52_struc{y}))/length(   apost_48_52_struc{y});
end



% CDF Real data

apre_9_11  = apre(apre<=11 & apre>= 9);
apre_48_52 = apre(apre<=52 & apre>= 48);
apost_9_11 = apost(apost<=11 & apost>= 9);
x_pre_50  = (1:length(apre_48_52))/length(apre_48_52);
x_pre   = (1:length(apre_9_11))/length(apre_9_11);

x_post  = (1:length(apost_9_11))/length(apost_9_11);

% 10B
figure;
plot(sort(apost_9_11),x_post,'-black','LineWidth',2)
hold on;
plot(sort(apost_9_11_struc{1}), x_post_struc{1} ,':','LineWidth',1.3,"Color", [0, 0.4470, 0.7410, 0.5])
plot(sort(apost_9_11_struc{2}), x_post_struc{2}  ,':','LineWidth',1.3,"Color", [0.8500, 0.3250, 0.0980, 0.5])
plot(sort(apost_9_11_struc{3}), x_post_struc{3}  ,':','LineWidth',1.3,"Color", [0.9290, 0.6940, 0.1250, 0.5])
plot(sort(apost_9_11_struc{4}), x_post_struc{4}  ,':','LineWidth',1.3,"Color", [0.4940, 0.1840, 0.5560, 0.5])
plot(sort(apost_9_11_struc{5}), x_post_struc{5}  ,':','LineWidth',1.3,"Color", [0.4660, 0.6740, 0.1880, 0.5])
plot(sort(apost_9_11_struc{6}), x_post_struc{6}  ,':','LineWidth',1.3,"Color", [0.3010, 0.7450, 0.9330])
plot(sort(apost_9_11_struc{7}), x_post_struc{7}  ,':','LineWidth',1.3,"Color", [0.6350, 0.0780, 0.184, 0.5])
plot([10 10],ylim,'-k')
hold off;
leg1 = legend('Post Dodd-Frank', '2011','2012','2013','2014','2015','2016','2017');
set(leg1,'Location','southeast');
xlabel('Assets')
fig = gcf;
fig.PaperPositionMode =  'manual'; %'auto'
fig.PaperPosition     = [1.3240    3.3281    5.8521    4.3438];
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
print(fig,'asset_distribution_cumulative_by_year_postDF_10B_designed','-dpdf','-fillpage')


%% Figures OA.2.: CDF for each year pre DF: asset_distribution_cumulative_by_year_preDF_10B_designed_until_2010q1.pdf

amin = 3;
amax = 250;
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.datequarter_num >= 2001 ); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = data(data.datequarter_num >= 2001 & data.datequarter_num <= 2010.25,:).assets;
apost = data(data.post == 2,:).assets;


data_pre = data(data.datequarter_num >= 2001 & data.datequarter_num <= 2010.25, :);
data_pre.year = floor(data_pre.datequarter_num);
data_pre = data_pre(data_pre.year >=2001 & data_pre.year <= 2010, :);

yearlist = sort(unique(data_pre.year));

apre_9_11_struc = {};
apost_9_11_struc = {};
apost_48_52_struc = {};
x_pre_struc = {};
x_pre_struct_50b = {};
x_post_struc = {};
x_post_struct_50b = {};


for y = 1:length(yearlist)
    apre_9_11_struc{y} = data_pre(data_pre.assets <=11 & data_pre.assets>= 9 & data_pre.year == yearlist(y),:).assets;
    apre_48_52_struc{y} = data_pre(data_pre.assets <=52 & data_pre.assets>= 48 & data_pre.year == yearlist(y),:).assets;

    yearlist(y)
    length(apre_9_11_struc{y})
    x_pre_struc{y}     = (1:length( apre_9_11_struc{y}))/length( apre_9_11_struc{y});
    x_pre_struct_50b{y} = (1:length(  apre_48_52_struc{y}))/length(   apre_48_52_struc{y});
end



% CDF Real data

apre_9_11  = apre(apre<=11 & apre>= 9);
apre_48_52 = apre(apre<=52 & apre>= 48);
apost_9_11 = apost(apost<=11 & apost>= 9);
x_pre_50  = (1:length(apre_48_52))/length(apre_48_52);
x_pre   = (1:length(apre_9_11))/length(apre_9_11);

x_post  = (1:length(apost_9_11))/length(apost_9_11);


% 10B
figure;
plot(sort(apre_9_11),x_pre,'-black','LineWidth',2)
hold on;
plot(sort(apre_9_11_struc{1}), x_pre_struc{1} ,':','LineWidth',1.3,"Color", [0, 0.4470, 0.7410, 0.5])
plot(sort(apre_9_11_struc{2}), x_pre_struc{2}  ,':','LineWidth',1.3,"Color", [0.8500, 0.3250, 0.0980, 0.5])
plot(sort(apre_9_11_struc{3}), x_pre_struc{3}  ,':','LineWidth',1.3,"Color", [0.9290, 0.6940, 0.1250, 0.5])
plot(sort(apre_9_11_struc{4}), x_pre_struc{4}  ,':','LineWidth',1.3,"Color", [0.4940, 0.1840, 0.5560, 0.5])
plot(sort(apre_9_11_struc{5}), x_pre_struc{5}  ,':','LineWidth',1.3,"Color", [0.4660, 0.6740, 0.1880, 0.5])
plot(sort(apre_9_11_struc{6}), x_pre_struc{6}  ,':','LineWidth',1.3,"Color", [0.3010, 0.7450, 0.9330])
plot(sort(apre_9_11_struc{7}), x_pre_struc{7}  ,':','LineWidth',1.3,"Color", [0.6350, 0.0780, 0.184, 0.5])
plot(sort(apre_9_11_struc{8}), x_pre_struc{8}  ,':','LineWidth',1.3,"Color", [0, 1.0000, 0.6667, 0.5])
plot(sort(apre_9_11_struc{9}), x_pre_struc{9}  ,':','LineWidth',1.3,"Color", [0.1804,    0.1804,    0.7882, 0.5])
plot(sort(apre_9_11_struc{10}), x_pre_struc{10}  ,':','LineWidth',1.3,"Color", [0.6118  ,  0.1647    ,0.1647, 0.5])
%plot(sort(apost_9_11_struc{8}), x_post_struc{8} ,'--k','LineWidth',2)
plot([10 10],ylim,'-k')
hold off;
leg1 = legend('Pre Dodd-Frank', '2001','2002','2003','2004','2005','2006','2007','2008','2009','2010 Q1');
set(leg1,'Location','southeast');
xlabel('Assets')
fig = gcf;
fig.PaperPositionMode =  'manual'; %'auto'
fig.PaperPosition     = [1.3240    3.3281    5.8521    4.3438];
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
print(fig,'asset_distribution_cumulative_by_year_preDF_10B_designed_until_2010q1_v2','-dpdf','-fillpage')

%% Figure 3: Illustrative Graphs (Theoretical effect of regulation on assets and profit)
k = 0.01/(10*0.03);
q_bar = 10;
q_up = exp(q_upperbar_proportional(k, q_bar ,[]));

figure;
x  = 7:0.001:16;
afun = @(q) (q<=q_bar) .* q + (q>q_bar) .* (q<=q_up) .* q_bar + (q>q_up) .* q;

hold on;
plot(x,afun(x) ,'-','LineWidth',2)
plot(x,x ,'--','LineWidth',2)
plot([q_bar,q_up;q_bar,q_up],[ylim;ylim]','--k')
% xticks([2,2.2,2.4,2.6,2.8])
% xticklabels({'log(7.4)','log(9)','log(11)','log(13.5)','log(16.4)'})
% yticks([2,2.2,2.4,2.6,2.8])
% yticklabels({'log(7.4)','log(9)','log(11)','log(13.5)','log(16.4)'})
xticks([])
yticks([])
ylabel('Post-regulation assets')
xlabel('Pre-regulation assets')
% legend('k = 0.01','Location','NorthWest')
% title('Assets')
grid off
text([q_bar+0.03 q_up+0.03], 14*[1 1], {'   $\underline{q}$', '   $\overline{q}$'},'Interpreter', 'latex')
text([7.4 7.5 q_bar q_up-0.2], [7.3 9.5 9.5 12], {') 45^\circ',' No effect', '            Bunching', '    Incur regulatory cost'})
hold off


fig = gcf;
fig.PaperPositionMode =  'manual'; %'auto'
fig.PaperPosition     = [1.3240    3.3281    5.8521    4.3438];
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
print(fig,'effect_of_regulation_on_assets','-dpdf','-fillpage')

R   = 0.05;
bet = 66;
pifun = @(R,bet,q,q_bar,q_up,k) (k>0) .* ((q<=q_bar) .* bet^-1 .* q + (q>q_bar) .* (q<=q_up) .* ( R - (bet^-1) .* (log(q_bar) - log(q) + bet.*R-1) ).*q_bar + (q>q_up) .* bet^-1 .* q .*(1-k)) + (k==0) .* bet^-1 .* q  ;
x  = 7:0.001:16;


figure;
hold on;
plot([q_bar:0.01:max(x)],(pifun(R,bet,[q_bar:0.01:max(x)],q_bar,max(x),0) - pifun(R,bet,[q_bar:0.01:max(x)],q_bar,max(x),k))./pifun(R,bet,[q_bar:0.01:max(x)],q_bar,max(x),0) ,'--','LineWidth',2)
plot([q_bar+0.001:0.01:max(x)],(pifun(R,bet,[q_bar+0.001:0.01:max(x)],q_bar,q_bar,0) - pifun(R,bet,[q_bar+0.001:0.01:max(x)],q_bar,q_bar,k))./pifun(R,bet,[q_bar+0.001:0.01:max(x)],q_bar,q_bar,0) ,':','LineWidth',2)
plot(x,(pifun(R,bet,x,q_bar,q_up,0) - pifun(R,bet,x,q_bar,q_up,k))./pifun(R,bet,x,q_bar,q_up,0) ,'-','LineWidth',2)
plot([q_bar,q_up;q_bar,q_up],[ylim;ylim]','--k')
ylabel('Profit loss')
xlabel('Pre-regulation assets')
leg1 = legend( 'Cost of bunching','Regulatory cost','Loss in profit');
set(leg1,'Location','northwest');
xticks([])
yticks([])
grid off
text([q_bar+0.04 q_up+0.04], 0.005*[1 1], {'   $\underline{q}$', '   $\overline{q}$'},'Interpreter', 'latex')
text([7.5 q_bar q_up-0.2], [20e-3 20e-3 20e-3], {'No effect', '            Bunching', '    Incur regulatory cost'})
hold off



fig = gcf;
fig.PaperPositionMode =  'manual'; %'auto'
fig.PaperPosition     = [1.3240    3.3281    5.8521    4.3438];
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
print(fig,'effect_of_regulation_on_profit','-dpdf','-fillpage')

%% Figure 5: Illustrative Graphs (Theoretical Asset distribution CDF )
q_bar = 10;
qmin = 0.001;
x = [8:0.2:9.8,9.9, 9.95, 9.99,9.99:0.001:10.001,10.002:0.1:16];
betaq = 1.15;

% CDF
figure;
subplot(2,2,1)
hold on
plot(x, F_distribution_power(x, q_bar, qmin, betaq, 0, exp(q_upperbar_proportional(0, q_bar ,[])))  ,'-','LineWidth',2)
plot(x, F_distribution_power(x, q_bar, qmin, betaq, 0, exp(q_upperbar_proportional(0.025, q_bar ,[]))),'--','LineWidth',2)
plot(x, F_distribution_power(x, q_bar, qmin, betaq, 0, exp(q_upperbar_proportional(0.05, q_bar ,[]))),':','LineWidth',2)
ylabel('Fraction of observation')
xlabel('Assets')
% legend('\tau = 0','\tau = 0.025','\tau = 0.050','Location','SouthEast')
title('\sigma = 0')
xlim([min(x) max(x)])
grid off
hold off

subplot(2,2,2)
hold on
plot(x, F_distribution_power(x, q_bar, qmin, betaq, 0.025, exp(q_upperbar_proportional(0, q_bar ,[])))  ,'-','LineWidth',2)
plot(x, F_distribution_power(x, q_bar, qmin, betaq, 0.025, exp(q_upperbar_proportional(0.025, q_bar ,[]))),'--','LineWidth',2)
plot(x, F_distribution_power(x, q_bar, qmin, betaq, 0.025, exp(q_upperbar_proportional(0.05, q_bar ,[]))),':','LineWidth',2)
ylabel('Fraction of observation')
xlabel('Assets')
% legend('\tau = 0','\tau = 0.025','\tau = 0.050','Location','SouthEast')
title('\sigma = 0.025') 
xlim([min(x) max(x)])
grid off
hold off

subplot(2,2,3)
hold on
plot(x, F_distribution_power(x, q_bar, qmin, betaq, 0.05, exp(q_upperbar_proportional(0, q_bar ,[])))  ,'-','LineWidth',2)
plot(x, F_distribution_power(x, q_bar, qmin, betaq, 0.05, exp(q_upperbar_proportional(0.025, q_bar ,[]))),'--','LineWidth',2)
plot(x, F_distribution_power(x, q_bar, qmin, betaq, 0.05, exp(q_upperbar_proportional(0.05, q_bar ,[]))),':','LineWidth',2)
ylabel('Fraction of observation')
xlabel('Assets')
% legend('\tau = 0','\tau = 0.025','\tau = 0.050','Location','SouthEast')
title('\sigma = 0.05')
xlim([min(x) max(x)])
grid off
hold off

subplot(2,2,4)
hold on
plot(x, F_distribution_power(x, q_bar, qmin, betaq, 0.1, exp(q_upperbar_proportional(0, q_bar ,[])))  ,'-','LineWidth',2)
plot(x, F_distribution_power(x, q_bar, qmin, betaq, 0.1, exp(q_upperbar_proportional(0.025, q_bar ,[]))),'--','LineWidth',2)
plot(x, F_distribution_power(x, q_bar, qmin, betaq, 0.1, exp(q_upperbar_proportional(0.05, q_bar ,[]))),':','LineWidth',2)
ylabel('Fraction of observation')
xlabel('Assets')
% legend('\tau = 0','\tau = 0.025','\tau = 0.050','Location','SouthEast')
title('\sigma = 0.10')
xlim([min(x) max(x)])
grid off
hold off

fig = gcf;
fig.PaperPositionMode =  'manual'; %'auto'
fig.PaperPosition     = [1.3240    3.3281    5.8521    4.3438];
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
print(fig,'theoretical_asset_cdf_vs_sigma_and_k_proportional_power_law','-dpdf','-fillpage')


%% Figure 6 a) Illustrative Graphs (Simulated versus actual data ) $10 Billion Threshold

rng(5); % 5
qmin   = 0.001;
betaq = 1.114;
sig   = 0.0426;
k     = 0.00406;
q_bar = 10;
% actual data
amin = 3;
amax = 25;
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = data(data.post == 0,:).assets;
apost = data(data.post == 2,:).assets;


% sim data
% Finv = @(c,betaq,qmin,F) (F * (1-betaq)/c + qmin ^(1-betaq)).^(1/(1-betaq));
Finv = @(c,betaq,qmin,F) ((1-F) .* ((betaq-1)/c)).^(1/(1-betaq));

c = (betaq-1)/qmin^(1-betaq);

NB = height(data_orig((data_orig.post==2 | data_orig.post==0) & data_orig.assets>=1,:));
% NB = 10000000
% NB = 100000; %round((length(apre)+length(apost))/2) ;
v     =  lognrnd(0,sig,NB,1);
astar = Finv(c,betaq,qmin,rand(NB,1));

apre_sim = astar.*v;
apre_sim = apre_sim(apre_sim>= 9 & apre_sim <= 11);

v     =  lognrnd(0,sig,NB,1);
astar = Finv(c,betaq,qmin,rand(NB,1));

astar_post = astar;
astar_post(astar>=q_bar & astar <= exp(q_upperbar_proportional(k, q_bar ,[]))) = q_bar;
apost_sim = astar_post.*v;
apost_sim = apost_sim(apost_sim>= 9 & apost_sim <= 11);

x_pre_sim   = (1:length(apre_sim))/length(apre_sim);
x_post_sim  = (1:length(apost_sim))/length(apost_sim);


% Actual
apre_9_11 = apre(apre<=11 & apre>= 9);
apost_9_11 = apost(apost<=11 & apost>= 9);
x_pre   = (1:length(apre_9_11))/length(apre_9_11);
x_post  = (1:length(apost_9_11))/length(apost_9_11);


% Both together
figure;
plot(sort(apre_sim),x_pre_sim ,'-k','Color',[0 0.4470 0.7410],'LineWidth',2)
hold on;
plot(sort(apre_9_11),x_pre ,'--','Color',[0 0.4470 0.7410],'LineWidth',2)
plot(sort(apost_sim),x_post_sim,'-','Color',[0.8500 0.3250 0.0980] ,'LineWidth',2)
plot(sort(apost_9_11),x_post ,'--','Color',[0.8500 0.3250 0.0980],'LineWidth',2)
%plot([9:0.1:11],[0:0.05:1],':','LineWidth',2)
plot([10 10],ylim,'-k')
% xlim([log(9), log(11)])
% xticklabels({'log(9)','log(9.5)','log(10)','log(10.5)','log(11)'})
hold off;
%title('Simulated Data')
legend('Model (Pre DF)', 'Data (Pre DF)','Model (Post DF)','Data (Post DF)','Location','SouthEast')
xlabel('Assets')
ylabel('Fraction of observations')
grid off
fig = gcf;
fig.PaperPositionMode = 'manual'; %'auto'
fig.PaperPosition     = [1.7417    3.5917    5.0167    3.8167];
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
print(fig,'asset_distribution_simulation_versus_data_actual_pre_post_10B','-dpdf','-fillpage')


%% Figure 6 b) Illustrative Graphs (Simulated versus actual data ) $50 Billion Threshold

rng(112);

qmin   = 0.001;
betaq = 1.083;
sig   = 0.0229;
k     = 0.00106;
q_bar = 50;
% actual data
amin = 25;
amax = 250;
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = data(data.post == 0,:).assets;
apost = data(data.post == 2,:).assets;

amin_gt = 48;
amax_gt = 52;

% sim data
% Finv = @(c,betaq,qmin,F) (F * (1-betaq)/c + qmin ^(1-betaq)).^(1/(1-betaq));
Finv = @(c,betaq,qmin,F) ((1-F) .* ((betaq-1)/c)).^(1/(1-betaq));

c = (betaq-1)/qmin^(1-betaq);

NB = height(data_orig((data_orig.post==2 | data_orig.post==0) & data_orig.assets>=1,:));
%NB = 10000000

% NB = 130000; %round((length(apre)+length(apost))/2) ;
v     =  lognrnd(0,sig,NB,1);
astar = Finv(c,betaq,qmin,rand(NB,1));

apre_sim = astar.*v;
apre_sim = apre_sim(apre_sim>= amin_gt & apre_sim <= amax_gt);

v     =  lognrnd(0,sig,NB,1);
astar = Finv(c,betaq,qmin,rand(NB,1));

astar_post = astar;
astar_post(astar>=q_bar & astar <= exp(q_upperbar_proportional(k, q_bar ,[]))) = q_bar;
apost_sim = astar_post.*v;
apost_sim = apost_sim(apost_sim>= amin_gt & apost_sim <= amax_gt);

x_pre_sim   = (1:length(apre_sim))/length(apre_sim);
x_post_sim  = (1:length(apost_sim))/length(apost_sim);

% Actual
apre_9_11 = apre(apre<=amax_gt & apre>= amin_gt);
apost_9_11 = apost(apost<=amax_gt & apost>= amin_gt);
x_pre   = (1:length(apre_9_11))/length(apre_9_11);
x_post  = (1:length(apost_9_11))/length(apost_9_11);


% Both together
figure;
plot(sort(apre_sim),x_pre_sim ,'-k','Color',[0 0.4470 0.7410],'LineWidth',2)
hold on;
plot(sort(apre_9_11),x_pre ,'--','Color',[0 0.4470 0.7410],'LineWidth',2)
plot(sort(apost_sim),x_post_sim,'-','Color',[0.8500 0.3250 0.0980] ,'LineWidth',2)
plot(sort(apost_9_11),x_post ,'--','Color',[0.8500 0.3250 0.0980],'LineWidth',2)
%plot([9:0.1:11],[0:0.05:1],':','LineWidth',2)
plot([50 50],ylim,'-k')
% xlim([log(9), log(11)])
% xticklabels({'log(9)','log(9.5)','log(10)','log(10.5)','log(11)'})
hold off;
%title('Simulated Data')
legend('Model (Pre DF)', 'Data (Pre DF)','Model (Post DF)','Data (Post DF)','Location','SouthEast')
xlabel('Assets')
ylabel('Fraction of observations')
grid off
fig = gcf;
fig.PaperPositionMode = 'manual'; %'auto'
fig.PaperPosition     = [1.7417    3.5917    5.0167    3.8167];
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
print(fig,'asset_distribution_simulation_versus_data_actual_pre_post_50B','-dpdf','-fillpage')



%% Table 2, Panel A: $10 Billion Threshold (Pre and Post Dodd-Frank):

% Starting with correct values: -------------------------------------------
amin = 3;
amax = 40;
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = []; %data(data.post == 0,:).assets;
apost = data(data.post == 2,:).assets;

EVENT_QUARTER = 2010.5;
START_QUARTER = min(data.datequarter_num);
data.t = zeros(height(data),1);
data(data.datequarter_num <EVENT_QUARTER,:).t = data(data.datequarter_num <EVENT_QUARTER,:).datequarter_num -  START_QUARTER;
data(data.datequarter_num >=EVENT_QUARTER,:).t = data(data.datequarter_num >=EVENT_QUARTER,:).datequarter_num -  EVENT_QUARTER;
data.IS_POST = data.post == 2;
nb_muz = length(unique(data.datequarter_num));

groupsummary(data,'datequarter_num','mean','assets')
groupsummary(data,'datequarter_num')
data2 = groupsummary(data,'datequarter_num');

data3 = join(data(:, {'datequarter_num','post','rssdid'}), data2);

data2 = sortrows(data,6);
% groupsummary(data,'datequarter_num','mean','assets')

q_bar = 10;
qmin   = 0.001;



% loglik_f_power_with_trend(data,EVENT_QUARTER, q_bar, qmin, betaq, sig, k, mu0, mu1)
% loglik  = @(b) loglik_f_power_with_trend(data, EVENT_QUARTER, q_bar, qmin, b(1), b(2), b(3), b(4), b(5));
% loglik  = @(b) loglik_f_power_with_varying_betas(data, EVENT_QUARTER, q_bar, qmin, b(1:nb_muz),b(nb_muz+1), b(nb_muz+2), b(nb_muz+3));
% ( apost,apre, q_bar, qmin, betaq, sig, k)
betaq = 1.6;
c     = 0.6;
q_up  = 11;
mu0   = 0;
mu1   = 0;


A       =  [];
b       = [];
Aeq     = [];
beq     = [];
nonlcon = [];
options = optimset('Diagnostics','on','Display','iter', ...
                    'TolX',1e-30,'TolFun',1e-30,'MaxIter',1000,'Algorithm','interior-point','MaxFunEvals',1000);

% with k as the parameter
loglik  = @(b) loglik_f_power(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;0.01];
xstart  = x0;
lb      = [1e-9;1e-9;0];
up      = [Inf;0.1;Inf ];
[theta,stderr,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% with q_up as the parameter
loglik_f  = @(apost,apre, q_bar, qmin, betaq, sig, q_up) [log(f_density_power(apre, q_bar, qmin, betaq, sig, q_bar));... 
    log(f_density_power(apost, q_bar, qmin, betaq, sig, q_up))];

loglik  = @(b) loglik_f(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;q_up];
xstart  = x0;
lb      = [1e-9;1e-9;q_bar];
up      = [Inf;0.1;Inf ];
[theta2,stderr2,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% Export table
FID = fopen(fullfile('parameters_estimates_baseline_MLE_10B_powerlaw.tex'), 'w');
fprintf(FID, '\\begin{tabular*}{\\textwidth}{l@{\\extracolsep{\\fill}}lll@{}}		\\hline   &  & Estimated value & S.E. \\\\ \\hline ');
fprintf(FID, concatenate_text(3,{'$\\beta$  & Exponent of the power law distribution    & ',theta(1),'     & [',stderr(1),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\sigma$ &  Measurement error volatility (in \\%%)  & ',theta(2)*100,' & [',stderr(2)*100,'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\exp(\\overline{q})$  &  Assets of marginal bank (\\$ Billion) & ',theta2(3),'     & [',stderr2(3),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{' $\\tau$  & Cost of regulation (\\%% of profit)      & ',theta(3) * 100,'& [',stderr(3) * 100,'] \\\\ '}));
fprintf(FID, ' \\hline \\end{tabular*} ');


%% Table 12 Panel A: $10 Billion Threshold (Pre and Post Dodd-Frank) (Robustness: Removing 1 year):
drop_n_years = 1;

% Starting with correct values: -------------------------------------------
amin = 3;
amax = 40;
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);

EVENT_QUARTER = 2010.5;
START_QUARTER = min(data.datequarter_num);

apre  = []; %data(data.post == 0,:).assets;
apost = data(data.post == 2 & data.datequarter_num >= EVENT_QUARTER + drop_n_years,:).assets;


data.t = zeros(height(data),1);
data(data.datequarter_num <EVENT_QUARTER,:).t = data(data.datequarter_num <EVENT_QUARTER,:).datequarter_num -  START_QUARTER;
data(data.datequarter_num >=EVENT_QUARTER,:).t = data(data.datequarter_num >=EVENT_QUARTER,:).datequarter_num -  EVENT_QUARTER;
data.IS_POST = data.post == 2;
nb_muz = length(unique(data.datequarter_num));

groupsummary(data,'datequarter_num','mean','assets')
groupsummary(data,'datequarter_num')
data2 = groupsummary(data,'datequarter_num');

data3 = join(data(:, {'datequarter_num','post','rssdid'}), data2);

data2 = sortrows(data,6);
% groupsummary(data,'datequarter_num','mean','assets')

q_bar = 10;
qmin   = 0.001;



% loglik_f_power_with_trend(data,EVENT_QUARTER, q_bar, qmin, betaq, sig, k, mu0, mu1)
% loglik  = @(b) loglik_f_power_with_trend(data, EVENT_QUARTER, q_bar, qmin, b(1), b(2), b(3), b(4), b(5));
% loglik  = @(b) loglik_f_power_with_varying_betas(data, EVENT_QUARTER, q_bar, qmin, b(1:nb_muz),b(nb_muz+1), b(nb_muz+2), b(nb_muz+3));
% ( apost,apre, q_bar, qmin, betaq, sig, k)
betaq = 1.6;
c     = 0.6;
q_up  = 11;
mu0   = 0;
mu1   = 0;


A       =  [];
b       = [];
Aeq     = [];
beq     = [];
nonlcon = [];
options = optimset('Diagnostics','on','Display','iter', ...
                    'TolX',1e-30,'TolFun',1e-30,'MaxIter',1000,'Algorithm','interior-point','MaxFunEvals',1000);

% with k as the parameter
loglik  = @(b) loglik_f_power(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;0.01];
xstart  = x0;
lb      = [1e-9;1e-9;0];
up      = [Inf;0.1;Inf ];
[theta,stderr,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% with q_up as the parameter
loglik_f  = @(apost,apre, q_bar, qmin, betaq, sig, q_up) [log(f_density_power(apre, q_bar, qmin, betaq, sig, q_bar));... 
    log(f_density_power(apost, q_bar, qmin, betaq, sig, q_up))];

loglik  = @(b) loglik_f(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;q_up];
xstart  = x0;
lb      = [1e-9;1e-9;q_bar];
up      = [Inf;0.1;Inf ];
[theta2,stderr2,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% Export table
FID = fopen(fullfile('parameters_estimates_baseline_MLE_10B_powerlaw_drop_1_year_postDF.tex'), 'w');
fprintf(FID, '\\begin{tabular*}{\\textwidth}{l@{\\extracolsep{\\fill}}lll@{}}		\\hline   &  & Estimated value & S.E. \\\\ \\hline ');
fprintf(FID, concatenate_text(3,{'$\\beta$  & Exponent of the power law distribution    & ',theta(1),'     & [',stderr(1),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\sigma$ &  Measurement error volatility (in \\%%)  & ',theta(2)*100,' & [',stderr(2)*100,'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\exp(\\overline{q})$  &  Assets of marginal bank (\\$ Billion) & ',theta2(3),'     & [',stderr2(3),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{' $\\tau$  & Cost of regulation (\\%% of profit)      & ',theta(3) * 100,'& [',stderr(3) * 100,'] \\\\ '}));
fprintf(FID, ' \\hline \\end{tabular*} ');



%% Table 12 Panel B: $10 Billion Threshold (Pre and Post Dodd-Frank) (Robustness: Removing 2 year):
drop_n_years = 2;

% Starting with correct values: -------------------------------------------
amin = 3;
amax = 40;
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);

EVENT_QUARTER = 2010.5;
START_QUARTER = min(data.datequarter_num);

apre  = []; %data(data.post == 0,:).assets;
apost = data(data.post == 2 & data.datequarter_num >= EVENT_QUARTER + drop_n_years,:).assets;


data.t = zeros(height(data),1);
data(data.datequarter_num <EVENT_QUARTER,:).t = data(data.datequarter_num <EVENT_QUARTER,:).datequarter_num -  START_QUARTER;
data(data.datequarter_num >=EVENT_QUARTER,:).t = data(data.datequarter_num >=EVENT_QUARTER,:).datequarter_num -  EVENT_QUARTER;
data.IS_POST = data.post == 2;
nb_muz = length(unique(data.datequarter_num));

groupsummary(data,'datequarter_num','mean','assets')
groupsummary(data,'datequarter_num')
data2 = groupsummary(data,'datequarter_num');

data3 = join(data(:, {'datequarter_num','post','rssdid'}), data2);

data2 = sortrows(data,6);
% groupsummary(data,'datequarter_num','mean','assets')

q_bar = 10;
qmin   = 0.001;



% loglik_f_power_with_trend(data,EVENT_QUARTER, q_bar, qmin, betaq, sig, k, mu0, mu1)
% loglik  = @(b) loglik_f_power_with_trend(data, EVENT_QUARTER, q_bar, qmin, b(1), b(2), b(3), b(4), b(5));
% loglik  = @(b) loglik_f_power_with_varying_betas(data, EVENT_QUARTER, q_bar, qmin, b(1:nb_muz),b(nb_muz+1), b(nb_muz+2), b(nb_muz+3));
% ( apost,apre, q_bar, qmin, betaq, sig, k)
betaq = 1.6;
c     = 0.6;
q_up  = 11;
mu0   = 0;
mu1   = 0;


A       =  [];
b       = [];
Aeq     = [];
beq     = [];
nonlcon = [];
options = optimset('Diagnostics','on','Display','iter', ...
                    'TolX',1e-30,'TolFun',1e-30,'MaxIter',1000,'Algorithm','interior-point','MaxFunEvals',1000);

% with k as the parameter
loglik  = @(b) loglik_f_power(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;0.01];
xstart  = x0;
lb      = [1e-9;1e-9;0];
up      = [Inf;0.1;Inf ];
[theta,stderr,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% with q_up as the parameter
loglik_f  = @(apost,apre, q_bar, qmin, betaq, sig, q_up) [log(f_density_power(apre, q_bar, qmin, betaq, sig, q_bar));... 
    log(f_density_power(apost, q_bar, qmin, betaq, sig, q_up))];

loglik  = @(b) loglik_f(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;q_up];
xstart  = x0;
lb      = [1e-9;1e-9;q_bar];
up      = [Inf;0.1;Inf ];
[theta2,stderr2,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% Export table
FID = fopen(fullfile('parameters_estimates_baseline_MLE_10B_powerlaw_drop_2_year_postDF.tex'), 'w');
fprintf(FID, '\\begin{tabular*}{\\textwidth}{l@{\\extracolsep{\\fill}}lll@{}}		\\hline   &  & Estimated value & S.E. \\\\ \\hline ');
fprintf(FID, concatenate_text(3,{'$\\beta$  & Exponent of the power law distribution    & ',theta(1),'     & [',stderr(1),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\sigma$ &  Measurement error volatility (in \\%%)  & ',theta(2)*100,' & [',stderr(2)*100,'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\exp(\\overline{q})$  &  Assets of marginal bank (\\$ Billion) & ',theta2(3),'     & [',stderr2(3),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{' $\\tau$  & Cost of regulation (\\%% of profit)      & ',theta(3) * 100,'& [',stderr(3) * 100,'] \\\\ '}));
fprintf(FID, ' \\hline \\end{tabular*} ');

%% Table 12 Panel C: $10 Billion Threshold (Pre and Post Dodd-Frank) (Robustness: Removing 3 year):
drop_n_years = 3;

% Starting with correct values: -------------------------------------------
amin = 3;
amax = 40;
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);

EVENT_QUARTER = 2010.5;
START_QUARTER = min(data.datequarter_num);

apre  = []; %data(data.post == 0,:).assets;
apost = data(data.post == 2 & data.datequarter_num >= EVENT_QUARTER + drop_n_years,:).assets;


data.t = zeros(height(data),1);
data(data.datequarter_num <EVENT_QUARTER,:).t = data(data.datequarter_num <EVENT_QUARTER,:).datequarter_num -  START_QUARTER;
data(data.datequarter_num >=EVENT_QUARTER,:).t = data(data.datequarter_num >=EVENT_QUARTER,:).datequarter_num -  EVENT_QUARTER;
data.IS_POST = data.post == 2;
nb_muz = length(unique(data.datequarter_num));

groupsummary(data,'datequarter_num','mean','assets')
groupsummary(data,'datequarter_num')
data2 = groupsummary(data,'datequarter_num');

data3 = join(data(:, {'datequarter_num','post','rssdid'}), data2);

data2 = sortrows(data,6);
% groupsummary(data,'datequarter_num','mean','assets')

q_bar = 10;
qmin   = 0.001;



% loglik_f_power_with_trend(data,EVENT_QUARTER, q_bar, qmin, betaq, sig, k, mu0, mu1)
% loglik  = @(b) loglik_f_power_with_trend(data, EVENT_QUARTER, q_bar, qmin, b(1), b(2), b(3), b(4), b(5));
% loglik  = @(b) loglik_f_power_with_varying_betas(data, EVENT_QUARTER, q_bar, qmin, b(1:nb_muz),b(nb_muz+1), b(nb_muz+2), b(nb_muz+3));
% ( apost,apre, q_bar, qmin, betaq, sig, k)
betaq = 1.6;
c     = 0.6;
q_up  = 11;
mu0   = 0;
mu1   = 0;


A       =  [];
b       = [];
Aeq     = [];
beq     = [];
nonlcon = [];
options = optimset('Diagnostics','on','Display','iter', ...
                    'TolX',1e-30,'TolFun',1e-30,'MaxIter',1000,'Algorithm','interior-point','MaxFunEvals',1000);

% with k as the parameter
loglik  = @(b) loglik_f_power(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;0.00458];
xstart  = x0;
lb      = [1e-9;1e-9;0];
up      = [Inf;0.1;Inf ];
[theta,stderr,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% with q_up as the parameter
loglik_f  = @(apost,apre, q_bar, qmin, betaq, sig, q_up) [log(f_density_power(apre, q_bar, qmin, betaq, sig, q_bar));... 
    log(f_density_power(apost, q_bar, qmin, betaq, sig, q_up))];

loglik  = @(b) loglik_f(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;q_up];
xstart  = x0;
lb      = [1e-9;1e-9;q_bar];
up      = [Inf;0.1;Inf ];
[theta2,stderr2,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% Export table
FID = fopen(fullfile('parameters_estimates_baseline_MLE_10B_powerlaw_drop_3_year_postDF.tex'), 'w');
fprintf(FID, '\\begin{tabular*}{\\textwidth}{l@{\\extracolsep{\\fill}}lll@{}}		\\hline   &  & Estimated value & S.E. \\\\ \\hline ');
fprintf(FID, concatenate_text(3,{'$\\beta$  & Exponent of the power law distribution    & ',theta(1),'     & [',stderr(1),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\sigma$ &  Measurement error volatility (in \\%%)  & ',theta(2)*100,' & [',stderr(2)*100,'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\exp(\\overline{q})$  &  Assets of marginal bank (\\$ Billion) & ',theta2(3),'     & [',stderr2(3),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{' $\\tau$  & Cost of regulation (\\%% of profit)      & ',theta(3) * 100,'& [',stderr(3) * 100,'] \\\\ '}));
fprintf(FID, ' \\hline \\end{tabular*} ');

%% Table OA.2.: $10 Billion Threshold (Pre and Post Dodd-Frank, 3-30 interval):

% Starting with correct values: -------------------------------------------
amin = 3;
amax = 30;
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = []; %data(data.post == 0,:).assets;
apost = data(data.post == 2,:).assets;

EVENT_QUARTER = 2010.5;
START_QUARTER = min(data.datequarter_num);
data.t = zeros(height(data),1);
data(data.datequarter_num <EVENT_QUARTER,:).t = data(data.datequarter_num <EVENT_QUARTER,:).datequarter_num -  START_QUARTER;
data(data.datequarter_num >=EVENT_QUARTER,:).t = data(data.datequarter_num >=EVENT_QUARTER,:).datequarter_num -  EVENT_QUARTER;
data.IS_POST = data.post == 2;
nb_muz = length(unique(data.datequarter_num));

groupsummary(data,'datequarter_num','mean','assets')
groupsummary(data,'datequarter_num')
data2 = groupsummary(data,'datequarter_num');

data3 = join(data(:, {'datequarter_num','post','rssdid'}), data2);

data2 = sortrows(data,6);
% groupsummary(data,'datequarter_num','mean','assets')

q_bar = 10;
qmin   = 0.001;



% loglik_f_power_with_trend(data,EVENT_QUARTER, q_bar, qmin, betaq, sig, k, mu0, mu1)
% loglik  = @(b) loglik_f_power_with_trend(data, EVENT_QUARTER, q_bar, qmin, b(1), b(2), b(3), b(4), b(5));
% loglik  = @(b) loglik_f_power_with_varying_betas(data, EVENT_QUARTER, q_bar, qmin, b(1:nb_muz),b(nb_muz+1), b(nb_muz+2), b(nb_muz+3));
% ( apost,apre, q_bar, qmin, betaq, sig, k)
betaq = 1.6;
c     = 0.6;
q_up  = 11;
mu0   = 0;
mu1   = 0;


A       =  [];
b       = [];
Aeq     = [];
beq     = [];
nonlcon = [];
options = optimset('Diagnostics','on','Display','iter', ...
                    'TolX',1e-30,'TolFun',1e-30,'MaxIter',1000,'Algorithm','interior-point','MaxFunEvals',1000);

% with k as the parameter
loglik  = @(b) loglik_f_power(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;0.01];
xstart  = x0;
lb      = [1e-9;1e-9;0];
up      = [Inf;0.1;Inf ];
[theta,stderr,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% with q_up as the parameter
loglik_f  = @(apost,apre, q_bar, qmin, betaq, sig, q_up) [log(f_density_power(apre, q_bar, qmin, betaq, sig, q_bar));... 
    log(f_density_power(apost, q_bar, qmin, betaq, sig, q_up))];

loglik  = @(b) loglik_f(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;q_up];
xstart  = x0;
lb      = [1e-9;1e-9;q_bar];
up      = [Inf;0.1;Inf ];
[theta2,stderr2,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% Export table
FID = fopen(fullfile('parameters_estimates_baseline_MLE_10B_powerlaw_3_30_interval.tex'), 'w');
fprintf(FID, '\\begin{tabular*}{\\textwidth}{l@{\\extracolsep{\\fill}}lll@{}}		\\hline   &  & Estimated value & S.E. \\\\ \\hline ');
fprintf(FID, concatenate_text(3,{'$\\beta$  & Exponent of the power law distribution    & ',theta(1),'     & [',stderr(1),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\sigma$ &  Measurement error volatility (in \\%%)  & ',theta(2)*100,' & [',stderr(2)*100,'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\exp(\\overline{q})$  &  Assets of marginal bank (\\$ Billion) & ',theta2(3),'     & [',stderr2(3),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{' $\\tau$  & Cost of regulation (\\%% of profit)      & ',theta(3) * 100,'& [',stderr(3) * 100,'] \\\\ '}));
fprintf(FID, ' \\hline \\end{tabular*} ');

%% Table 2, Panel B: $50 Billion Threshold (Pre and Post Dodd-Frank):

% Starting with correct values: -------------------------------------------
amin = 40;
amax = Inf;
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = [];% data(data.post == 0,:).assets;
apost = data(data.post == 2,:).assets;

EVENT_QUARTER = 2010.5;
START_QUARTER = min(data.datequarter_num);
data.t = zeros(height(data),1);
data(data.datequarter_num <EVENT_QUARTER,:).t = data(data.datequarter_num <EVENT_QUARTER,:).datequarter_num -  START_QUARTER;
data(data.datequarter_num >=EVENT_QUARTER,:).t = data(data.datequarter_num >=EVENT_QUARTER,:).datequarter_num -  EVENT_QUARTER;
data.IS_POST = data.post == 2;
nb_muz = length(unique(data.datequarter_num));

groupsummary(data,'datequarter_num','mean','assets')
groupsummary(data,'datequarter_num')
data2 = groupsummary(data,'datequarter_num');

data3 = join(data(:, {'datequarter_num','post','rssdid'}), data2);

data2 = sortrows(data,6);
% groupsummary(data,'datequarter_num','mean','assets')

q_bar = 50;
qmin   = 0.001;

% loglik_f_power_with_trend(data,EVENT_QUARTER, q_bar, qmin, betaq, sig, k, mu0, mu1)
% loglik  = @(b) loglik_f_power_with_trend(data, EVENT_QUARTER, q_bar, qmin, b(1), b(2), b(3), b(4), b(5));
% loglik  = @(b) loglik_f_power_with_varying_betas(data, EVENT_QUARTER, q_bar, qmin, b(1:nb_muz),b(nb_muz+1), b(nb_muz+2), b(nb_muz+3));
% ( apost,apre, q_bar, qmin, betaq, sig, k)
betaq = 1.6;
c     = 0.6;
q_up  = 11;
mu0   = 0;
mu1   = 0;


A       =  [];
b       = [];
Aeq     = [];
beq     = [];
nonlcon = [];
options = optimset('Diagnostics','on','Display','iter', ...
                    'TolX',1e-30,'TolFun',1e-30,'MaxIter',1000,'Algorithm','interior-point','MaxFunEvals',1000);

sig_up   = 0.05;
sig_down = 1e-9;
k_down    = 0;

% with k as the parameter
loglik  = @(b) loglik_f_power(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;0.01];
xstart  = x0;
lb      = [1e-9;sig_down;k_down];
up      = [Inf ;sig_up   ;Inf ];
[theta,stderr,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% with q_up as the parameter
loglik_f  = @(apost,apre, q_bar, qmin, betaq, sig, q_up) [log(f_density_power(apre, q_bar, qmin, betaq, sig, q_bar));... 
    log(f_density_power(apost, q_bar, qmin, betaq, sig, q_up))];

loglik  = @(b) loglik_f(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.04;q_up];
xstart  = x0;
lb      = [1e-9 ;sig_down ;q_bar];
up      = [Inf  ;sig_up   ;Inf ];
[theta2,stderr2,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% Export table
FID = fopen(fullfile('parameters_estimates_baseline_MLE_50B_powerlaw.tex'), 'w');
fprintf(FID, '\\begin{tabular*}{\\textwidth}{l@{\\extracolsep{\\fill}}lll@{}}		\\hline   &  & Estimated value & S.E. \\\\ \\hline ');
fprintf(FID, concatenate_text(3,{'$\\beta$  & Exponent of the power law distribution    & ',theta(1),'     & [',stderr(1),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\sigma$ &  Measurement error volatility (in \\%%)  & ',theta(2)*100,' & [',stderr(2)*100,'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\exp(\\overline{q})$  &  Assets of marginal bank (\\$ Billion) & ',theta2(3),'     & [',stderr2(3),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{' $\\tau$  & Cost of regulation (\\%% of profit)      & ',theta(3) * 100,'& [',stderr(3) * 100,'] \\\\ '}));
fprintf(FID, ' \\hline \\end{tabular*} ');

%% Table OA 3, Panel A) $10 Billion Threshold (Pre and Post Dodd-Frank) (Using pre regulation data):

% Starting with correct values: -------------------------------------------
amin = 3;
amax = 40;
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = data(data.post == 0,:).assets;
apost = data(data.post == 2,:).assets;

EVENT_QUARTER = 2010.5;
START_QUARTER = min(data.datequarter_num);
data.t = zeros(height(data),1);
data(data.datequarter_num <EVENT_QUARTER,:).t = data(data.datequarter_num <EVENT_QUARTER,:).datequarter_num -  START_QUARTER;
data(data.datequarter_num >=EVENT_QUARTER,:).t = data(data.datequarter_num >=EVENT_QUARTER,:).datequarter_num -  EVENT_QUARTER;
data.IS_POST = data.post == 2;
nb_muz = length(unique(data.datequarter_num));

groupsummary(data,'datequarter_num','mean','assets')
groupsummary(data,'datequarter_num')
data2 = groupsummary(data,'datequarter_num');

data3 = join(data(:, {'datequarter_num','post','rssdid'}), data2);

data2 = sortrows(data,6);
% groupsummary(data,'datequarter_num','mean','assets')

q_bar = 10;
qmin   = 0.001;



% loglik_f_power_with_trend(data,EVENT_QUARTER, q_bar, qmin, betaq, sig, k, mu0, mu1)
% loglik  = @(b) loglik_f_power_with_trend(data, EVENT_QUARTER, q_bar, qmin, b(1), b(2), b(3), b(4), b(5));
% loglik  = @(b) loglik_f_power_with_varying_betas(data, EVENT_QUARTER, q_bar, qmin, b(1:nb_muz),b(nb_muz+1), b(nb_muz+2), b(nb_muz+3));
% ( apost,apre, q_bar, qmin, betaq, sig, k)
betaq = 1.6;
c     = 0.6;
q_up  = 11;
mu0   = 0;
mu1   = 0;


A       =  [];
b       = [];
Aeq     = [];
beq     = [];
nonlcon = [];
options = optimset('Diagnostics','on','Display','iter', ...
                    'TolX',1e-30,'TolFun',1e-30,'MaxIter',1000,'Algorithm','interior-point','MaxFunEvals',1000);

% with k as the parameter
loglik  = @(b) loglik_f_power(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;0.01];
xstart  = x0;
lb      = [1e-9;1e-9;0];
up      = [Inf;0.1;Inf ];
[theta,stderr,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% with q_up as the parameter
loglik_f  = @(apost,apre, q_bar, qmin, betaq, sig, q_up) [log(f_density_power(apre, q_bar, qmin, betaq, sig, q_bar));... 
    log(f_density_power(apost, q_bar, qmin, betaq, sig, q_up))];

loglik  = @(b) loglik_f(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;q_up];
xstart  = x0;
lb      = [1e-9;1e-9;q_bar];
up      = [Inf;0.1;Inf ];
[theta2,stderr2,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% Export table
FID = fopen(fullfile('parameters_estimates_baseline_MLE_10B_powerlaw_using_pre_regul.tex'), 'w');
fprintf(FID, '\\begin{tabular*}{\\textwidth}{l@{\\extracolsep{\\fill}}lll@{}}		\\hline   &  & Estimated value & S.E. \\\\ \\hline ');
fprintf(FID, concatenate_text(3,{'$\\beta$  & Exponent of the power law distribution    & ',theta(1),'     & [',stderr(1),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\sigma$ &  Measurement error volatility (in \\%%)  & ',theta(2)*100,' & [',stderr(2)*100,'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\exp(\\overline{q})$  &  Assets of marginal bank (\\$ Billion) & ',theta2(3),'     & [',stderr2(3),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{' $\\tau$  & Cost of regulation (\\%% of profit)      & ',theta(3) * 100,'& [',stderr(3) * 100,'] \\\\ '}));
fprintf(FID, ' \\hline \\end{tabular*} ');


%% Table OA 3, Panel B) $50 Billion Threshold (Pre and Post Dodd-Frank) (using pre DF regulation data):

% Starting with correct values: -------------------------------------------
amin = 40;
amax = Inf;
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = data(data.post == 0,:).assets;
apost = data(data.post == 2,:).assets;

EVENT_QUARTER = 2010.5;
START_QUARTER = min(data.datequarter_num);
data.t = zeros(height(data),1);
data(data.datequarter_num <EVENT_QUARTER,:).t = data(data.datequarter_num <EVENT_QUARTER,:).datequarter_num -  START_QUARTER;
data(data.datequarter_num >=EVENT_QUARTER,:).t = data(data.datequarter_num >=EVENT_QUARTER,:).datequarter_num -  EVENT_QUARTER;
data.IS_POST = data.post == 2;
nb_muz = length(unique(data.datequarter_num));

groupsummary(data,'datequarter_num','mean','assets')
groupsummary(data,'datequarter_num')
data2 = groupsummary(data,'datequarter_num');

data3 = join(data(:, {'datequarter_num','post','rssdid'}), data2);

data2 = sortrows(data,6);
% groupsummary(data,'datequarter_num','mean','assets')

q_bar = 50;
qmin   = 0.001;

% loglik_f_power_with_trend(data,EVENT_QUARTER, q_bar, qmin, betaq, sig, k, mu0, mu1)
% loglik  = @(b) loglik_f_power_with_trend(data, EVENT_QUARTER, q_bar, qmin, b(1), b(2), b(3), b(4), b(5));
% loglik  = @(b) loglik_f_power_with_varying_betas(data, EVENT_QUARTER, q_bar, qmin, b(1:nb_muz),b(nb_muz+1), b(nb_muz+2), b(nb_muz+3));
% ( apost,apre, q_bar, qmin, betaq, sig, k)
betaq = 1.6;
c     = 0.6;
q_up  = 11;
mu0   = 0;
mu1   = 0;


A       =  [];
b       = [];
Aeq     = [];
beq     = [];
nonlcon = [];
options = optimset('Diagnostics','on','Display','iter', ...
                    'TolX',1e-30,'TolFun',1e-30,'MaxIter',1000,'Algorithm','interior-point','MaxFunEvals',1000);

sig_up   = 0.05;
sig_down = 1e-9;
k_down    = 0;

% with k as the parameter
loglik  = @(b) loglik_f_power(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;0.01];
xstart  = x0;
lb      = [1e-9;sig_down;k_down];
up      = [Inf ;sig_up   ;Inf ];
[theta,stderr,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% with q_up as the parameter
loglik_f  = @(apost,apre, q_bar, qmin, betaq, sig, q_up) [log(f_density_power(apre, q_bar, qmin, betaq, sig, q_bar));... 
    log(f_density_power(apost, q_bar, qmin, betaq, sig, q_up))];

loglik  = @(b) loglik_f(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.04;q_up];
xstart  = x0;
lb      = [1e-9 ;sig_down ;q_bar];
up      = [Inf  ;sig_up   ;Inf ];
[theta2,stderr2,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% Export table
FID = fopen(fullfile('parameters_estimates_baseline_MLE_50B_powerlaw_using_pre_regul.tex'), 'w');
fprintf(FID, '\\begin{tabular*}{\\textwidth}{l@{\\extracolsep{\\fill}}lll@{}}		\\hline   &  & Estimated value & S.E. \\\\ \\hline ');
fprintf(FID, concatenate_text(3,{'$\\beta$  & Exponent of the power law distribution    & ',theta(1),'     & [',stderr(1),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\sigma$ &  Measurement error volatility (in \\%%)  & ',theta(2)*100,' & [',stderr(2)*100,'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\exp(\\overline{q})$  &  Assets of marginal bank (\\$ Billion) & ',theta2(3),'     & [',stderr2(3),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{' $\\tau$  & Cost of regulation (\\%% of profit)      & ',theta(3) * 100,'& [',stderr(3) * 100,'] \\\\ '}));
fprintf(FID, ' \\hline \\end{tabular*} ');





%% Table OA.2, Panel B) $50 Billion Threshold (Pre and Post Dodd-Frank, 30 - Inf interval):

% Starting with correct values: -------------------------------------------
amin = 30;
amax = Inf;
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = []; %data(data.post == 0,:).assets;
apost = data(data.post == 2,:).assets;

EVENT_QUARTER = 2010.5;
START_QUARTER = min(data.datequarter_num);
data.t = zeros(height(data),1);
data(data.datequarter_num <EVENT_QUARTER,:).t = data(data.datequarter_num <EVENT_QUARTER,:).datequarter_num -  START_QUARTER;
data(data.datequarter_num >=EVENT_QUARTER,:).t = data(data.datequarter_num >=EVENT_QUARTER,:).datequarter_num -  EVENT_QUARTER;
data.IS_POST = data.post == 2;
nb_muz = length(unique(data.datequarter_num));

groupsummary(data,'datequarter_num','mean','assets')
groupsummary(data,'datequarter_num')
data2 = groupsummary(data,'datequarter_num');

data3 = join(data(:, {'datequarter_num','post','rssdid'}), data2);

data2 = sortrows(data,6);
% groupsummary(data,'datequarter_num','mean','assets')

q_bar = 50;
qmin   = 0.001;

% loglik_f_power_with_trend(data,EVENT_QUARTER, q_bar, qmin, betaq, sig, k, mu0, mu1)
% loglik  = @(b) loglik_f_power_with_trend(data, EVENT_QUARTER, q_bar, qmin, b(1), b(2), b(3), b(4), b(5));
% loglik  = @(b) loglik_f_power_with_varying_betas(data, EVENT_QUARTER, q_bar, qmin, b(1:nb_muz),b(nb_muz+1), b(nb_muz+2), b(nb_muz+3));
% ( apost,apre, q_bar, qmin, betaq, sig, k)
betaq = 1.6;
c     = 0.6;
q_up  = 11;
mu0   = 0;
mu1   = 0;


A       =  [];
b       = [];
Aeq     = [];
beq     = [];
nonlcon = [];
options = optimset('Diagnostics','on','Display','iter', ...
                    'TolX',1e-30,'TolFun',1e-30,'MaxIter',1000,'Algorithm','interior-point','MaxFunEvals',1000);

sig_up   = 0.05;
sig_down = 1e-9;
k_down    = 0;

% with k as the parameter
loglik  = @(b) loglik_f_power(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;0.01];
xstart  = x0;
lb      = [1e-9;sig_down;k_down];
up      = [Inf ;sig_up   ;Inf ];
[theta,stderr,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% with q_up as the parameter
loglik_f  = @(apost,apre, q_bar, qmin, betaq, sig, q_up) [log(f_density_power(apre, q_bar, qmin, betaq, sig, q_bar));... 
    log(f_density_power(apost, q_bar, qmin, betaq, sig, q_up))];

loglik  = @(b) loglik_f(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.04;q_up];
xstart  = x0;
lb      = [1e-9 ;sig_down ;q_bar];
up      = [Inf  ;sig_up   ;Inf ];
[theta2,stderr2,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% Export table
FID = fopen(fullfile('parameters_estimates_baseline_MLE_50B_powerlaw_30_inf_interval.tex'), 'w');
fprintf(FID, '\\begin{tabular*}{\\textwidth}{l@{\\extracolsep{\\fill}}lll@{}}		\\hline   &  & Estimated value & S.E. \\\\ \\hline ');
fprintf(FID, concatenate_text(3,{'$\\beta$  & Exponent of the power law distribution    & ',theta(1),'     & [',stderr(1),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\sigma$ &  Measurement error volatility (in \\%%)  & ',theta(2)*100,' & [',stderr(2)*100,'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\exp(\\overline{q})$  &  Assets of marginal bank (\\$ Billion) & ',theta2(3),'     & [',stderr2(3),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{' $\\tau$  & Cost of regulation (\\%% of profit)      & ',theta(3) * 100,'& [',stderr(3) * 100,'] \\\\ '}));
fprintf(FID, ' \\hline \\end{tabular*} ');



%% Table 9, Panel A) $10 Billion Threshold (Pre Dodd-Frank versus Post Regulatory Relief):

% Starting with correct values: -------------------------------------------
sigparam = 4.260/100;
betaq_param = 1.114;
amin = 3;
amax = 40;
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2 | data_orig.post == 3); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = []; % data(data.post == 0,:).assets;
apost = data(data.post == 3,:).assets;

EVENT_QUARTER = 2010.5;
START_QUARTER = min(data.datequarter_num);
data.t = zeros(height(data),1);
data(data.datequarter_num <EVENT_QUARTER,:).t = data(data.datequarter_num <EVENT_QUARTER,:).datequarter_num -  START_QUARTER;
data(data.datequarter_num >=EVENT_QUARTER,:).t = data(data.datequarter_num >=EVENT_QUARTER,:).datequarter_num -  EVENT_QUARTER;
data.IS_POST = data.post == 2;
nb_muz = length(unique(data.datequarter_num));

groupsummary(data,'datequarter_num','mean','assets')
groupsummary(data,'datequarter_num')
data2 = groupsummary(data,'datequarter_num');

data3 = join(data(:, {'datequarter_num','post','rssdid'}), data2);

data2 = sortrows(data,6);
% groupsummary(data,'datequarter_num','mean','assets')

q_bar = 10;
qmin   = 0.001;

% loglik_f_power_with_trend(data,EVENT_QUARTER, q_bar, qmin, betaq, sig, k, mu0, mu1)
% loglik  = @(b) loglik_f_power_with_trend(data, EVENT_QUARTER, q_bar, qmin, b(1), b(2), b(3), b(4), b(5));
% loglik  = @(b) loglik_f_power_with_varying_betas(data, EVENT_QUARTER, q_bar, qmin, b(1:nb_muz),b(nb_muz+1), b(nb_muz+2), b(nb_muz+3));
% ( apost,apre, q_bar, qmin, betaq, sig, k)
betaq = 1.6;
c     = 0.6;
q_up  = 11;
mu0   = 0;
mu1   = 0;


A       =  [];
b       = [];

nonlcon = [];
options = optimset('Diagnostics','on','Display','iter', ...
                    'TolX',1e-30,'TolFun',1e-30,'MaxIter',1000,'Algorithm','interior-point','MaxFunEvals',1000);

               
% with k as the parameter
loglik  = @(b) loglik_f_power(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;0.01];
xstart  = x0;
lb      = [1e-9;1e-10;0];
up      = [Inf;0.1;Inf ];
Aeq     = [1 0 0;0 1 0];
beq     = [betaq_param; sigparam];
[theta,stderr,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% with q_up as the parameter
loglik_f  = @(apost,apre, q_bar, qmin, betaq, sig, q_up) [log(f_density_power(apre, q_bar, qmin, betaq, sig, q_bar));... 
    log(f_density_power(apost, q_bar, qmin, betaq, sig, q_up))];

loglik  = @(b) loglik_f(apost, apre, q_bar, qmin, b(1),b(2), b(3));
x0      = [betaq;0.01;q_up];
xstart  = x0;
lb      = [1e-9;1e-10;q_bar];
up      = [Inf;0.1;Inf ];

[theta2,stderr2,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% Export table
FID = fopen(fullfile('parameters_estimates_baseline_MLE_10B_powerlaw_including_reg_relief.tex'), 'w');
fprintf(FID, '\\begin{tabular*}{\\textwidth}{l@{\\extracolsep{\\fill}}lll@{}}		\\hline   &  & Estimated value & S.E. \\\\ \\hline ');
%fprintf(FID, concatenate_text(3,{'$\\beta$  & Exponent of the power law distribution    & ',theta(1),'     & [','-','] \\\\ '}));
%fprintf(FID, concatenate_text(3,{'$\\sigma$ &  Measurement error volatility (in \\%%)  & ',theta(2)*100,' & [','-','] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\exp(\\overline{q})$  &  Assets of marginal bank (\\$ Billion) & ',theta2(3),'     & [',stderr2(3),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{' $\\tau$  & Cost of regulation (\\%% of profit)      & ',theta(3) * 100,'& [',stderr(3) * 100,'] \\\\ '}));
fprintf(FID, ' \\hline \\end{tabular*} ');



%% Table 9, Panel B) $50 Billion Threshold (Pre and Post Regulatory Relief):


% Starting with correct values: -------------------------------------------
sigparam = 2.29/100;
betaq_param = 1.090;
amin = 40;
amax = Inf;
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2 | data_orig.post == 3); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = data(data.post == 0,:).assets;
apost = data(data.post == 3,:).assets;

EVENT_QUARTER = 2010.5;
START_QUARTER = min(data.datequarter_num);
data.t = zeros(height(data),1);
data(data.datequarter_num <EVENT_QUARTER,:).t = data(data.datequarter_num <EVENT_QUARTER,:).datequarter_num -  START_QUARTER;
data(data.datequarter_num >=EVENT_QUARTER,:).t = data(data.datequarter_num >=EVENT_QUARTER,:).datequarter_num -  EVENT_QUARTER;
data.IS_POST = data.post == 2;
nb_muz = length(unique(data.datequarter_num));

groupsummary(data,'datequarter_num','mean','assets')
groupsummary(data,'datequarter_num')
data2 = groupsummary(data,'datequarter_num');

data3 = join(data(:, {'datequarter_num','post','rssdid'}), data2);

data2 = sortrows(data,6);
% groupsummary(data,'datequarter_num','mean','assets')

q_bar = 50;
qmin   = 0.001;

% loglik_f_power_with_trend(data,EVENT_QUARTER, q_bar, qmin, betaq, sig, k, mu0, mu1)
% loglik  = @(b) loglik_f_power_with_trend(data, EVENT_QUARTER, q_bar, qmin, b(1), b(2), b(3), b(4), b(5));
% loglik  = @(b) loglik_f_power_with_varying_betas(data, EVENT_QUARTER, q_bar, qmin, b(1:nb_muz),b(nb_muz+1), b(nb_muz+2), b(nb_muz+3));

% ( apost,apre, q_bar, qmin, betaq, sig, k)
betaq = 1.4;
c     = 0.6;
q_up  = 51;
mu0   = 0;
mu1   = 0;


A       =  [];
b       = [];
Aeq     = [];
beq     = [];
nonlcon = [];
options = optimset('Diagnostics','on','Display','iter', ...
                    'TolX',1e-30,'TolFun',1e-30,'MaxIter',1000,'Algorithm','interior-point','MaxFunEvals',1000);

% with k as the parameter
loglik  = @(b) loglik_f_power(apost, apre, q_bar, qmin, betaq_param,sigparam, b);
x0      = [0.001];
xstart  = x0;
lb      = [0];
up      = [Inf ];
[theta,stderr,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% with q_up as the parameter
loglik_f  = @(apost,apre, q_bar, qmin, betaq, sig, q_up) [log(f_density_power(apre, q_bar, qmin, betaq, sig, q_bar));... 
    log(f_density_power(apost, q_bar, qmin, betaq, sig, q_up))];

loglik  = @(b) loglik_f(apost, apre, q_bar, qmin, betaq_param,sigparam,  b);
x0      = [50]; % -28327.9333
xstart  = x0;
lb      = [q_bar];
up      = [Inf ];
[theta2,stderr2,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% Export table
FID = fopen(fullfile('parameters_estimates_baseline_MLE_50B_powerlaw_including_reg_relief.tex'), 'w');
fprintf(FID, '\\begin{tabular*}{\\textwidth}{l@{\\extracolsep{\\fill}}lll@{}}		\\hline   &  & Estimated value & S.E. \\\\ \\hline ');
%fprintf(FID, concatenate_text(3,{'$\\beta$  & Exponent of the power law distribution    & ',betaq_param,'     & [','-','] \\\\ '}));
%fprintf(FID, concatenate_text(3,{'$\\sigma$ &  Measurement error volatility (in \\%%)  & ',sigparam*100,' & [','-','] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\exp(\\overline{q})$  &  Assets of marginal bank (\\$ Billion) & ',theta2(1),'     & [',stderr2(1),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{' $\\tau$  & Cost of regulation (\\%% of profit)      & ',theta(1) * 100,'& [',stderr(1) * 100,'] \\\\ '}));
fprintf(FID, ' \\hline \\end{tabular*} ');



%% Table 11, Panel A): Placebo test Using Only Pre-DF data ($10B)


% Starting with correct values: -------------------------------------------
sigparam = 4.26/100;
betaq_param = 1.114;
amin = 3;
amax = 40;
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
% apre  = data(data.post == 0,:).assets;
% apost = data(data.post == 2,:).assets;

apost = data(data.post == 0,:).assets;
apre  = [];

EVENT_QUARTER = 2010.5;
START_QUARTER = min(data.datequarter_num);
data.t = zeros(height(data),1);
data(data.datequarter_num <EVENT_QUARTER,:).t = data(data.datequarter_num <EVENT_QUARTER,:).datequarter_num -  START_QUARTER;
data(data.datequarter_num >=EVENT_QUARTER,:).t = data(data.datequarter_num >=EVENT_QUARTER,:).datequarter_num -  EVENT_QUARTER;
data.IS_POST = data.post == 2;
nb_muz = length(unique(data.datequarter_num));

groupsummary(data,'datequarter_num','mean','assets')
groupsummary(data,'datequarter_num')
data2 = groupsummary(data,'datequarter_num');

data3 = join(data(:, {'datequarter_num','post','rssdid'}), data2);

data2 = sortrows(data,6);
% groupsummary(data,'datequarter_num','mean','assets')

q_bar = 10;
qmin   = 0.001;



% loglik_f_power_with_trend(data,EVENT_QUARTER, q_bar, qmin, betaq, sig, k, mu0, mu1)
% loglik  = @(b) loglik_f_power_with_trend(data, EVENT_QUARTER, q_bar, qmin, b(1), b(2), b(3), b(4), b(5));
% loglik  = @(b) loglik_f_power_with_varying_betas(data, EVENT_QUARTER, q_bar, qmin, b(1:nb_muz),b(nb_muz+1), b(nb_muz+2), b(nb_muz+3));
% ( apost,apre, q_bar, qmin, betaq, sig, k)
betaq = 1.6;
c     = 0.6;
q_up  = 11;
mu0   = 0;
mu1   = 0;


A       =  [];
b       = [];
Aeq     = [];
beq     = [];
nonlcon = [];
options = optimset('Diagnostics','on','Display','iter', ...
                    'TolX',1e-30,'TolFun',1e-30,'MaxIter',1000,'Algorithm','interior-point','MaxFunEvals',1000);

% with k as the parameter
loglik  = @(b) loglik_f_power(apost, apre, q_bar, qmin, b(1),sigparam, b(2));
x0      = [betaq;0.01];
xstart  = x0;
lb      = [1e-9;0];
up      = [Inf;Inf ];
[theta,stderr,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% with q_up as the parameter
loglik_f  = @(apost,apre, q_bar, qmin, betaq, sig, q_up) [log(f_density_power(apre, q_bar, qmin, betaq, sig, q_bar));... 
    log(f_density_power(apost, q_bar, qmin, betaq, sig, q_up))];

loglik  = @(b) loglik_f(apost, apre, q_bar, qmin, b(1),sigparam, b(2));
x0      = [betaq;q_up];
xstart  = x0;
lb      = [1e-9;q_bar];
up      = [Inf;Inf ];
[theta2,stderr2,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);



% 
% % Export table
% FID = fopen(fullfile(strcat(Tables,'/parameters_estimates_baseline_MLE_10B_powerlaw_preDF_placebo.tex')), 'w');
% fprintf(FID, '\\begin{tabular*}{\\textwidth}{l@{\\extracolsep{\\fill}}lll@{}}		\\hline   &  & \\textbf{Estimated value} & \\textbf{S.E.} \\\\ \\hline ');
% fprintf(FID, concatenate_text(3,{'$\\beta$  & Exponent of the power law distribution    & ',1.109,'     & [','-','] \\\\ '}));
% fprintf(FID, concatenate_text(3,{'$\\sigma$ &  Measurement error volatility (in \\%%)  & ',4.15,' & [','-','] \\\\ '}));
% fprintf(FID, concatenate_text(3,{'$\\exp(\\overline{q})$  &  Assets of marginal bank (\\$ Billion) & ',theta2(1),'     & [',stderr2(1),'] \\\\ '}));
% fprintf(FID, concatenate_text(3,{' $\\tau$  & Cost of regulation (\\%% of profit)      & ',theta(1) * 100,'& [',stderr(1) * 100,'] \\\\ '}));
% fprintf(FID, ' \\hline \\end{tabular*} ');



% Export table
FID = fopen(fullfile('parameters_estimates_baseline_MLE_10B_powerlaw_preDF_placebo.tex'), 'w');
fprintf(FID, '\\begin{tabular*}{\\textwidth}{l@{\\extracolsep{\\fill}}lll@{}}		\\hline   &  & Estimated value & S.E. \\\\ \\hline ');
fprintf(FID, concatenate_text(3,{'$\\beta$  & Exponent of the power law distribution    & ',theta(1),'     & [',stderr(1),'] \\\\ '}));
%fprintf(FID, concatenate_text(3,{'$\\sigma$ &  Measurement error volatility (in \\%%)  & ',sigparam*100,' & [','-','] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\exp(\\overline{q})$  &  Assets of marginal bank (\\$ Billion) & ',theta2(2),'     & [',stderr2(2),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{' $\\tau$  & Cost of regulation (\\%% of profit)      & ',theta(2) * 100,'& [',stderr(2) * 100,'] \\\\ '}));
fprintf(FID, ' \\hline \\end{tabular*} ');


%% Table 11, Panel B) Placebo test Using Only Pre-DF data ($50B)


% Starting with correct values: -------------------------------------------
sigparam = 2.29/100;
betaq_param = 1.09;
amin = 40;
amax = Inf;
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
% apre  = data(data.post == 0,:).assets;
% apost = data(data.post == 2,:).assets;

apost = data(data.post == 0,:).assets;
apre  = [];

EVENT_QUARTER = 2010.5;
START_QUARTER = min(data.datequarter_num);
data.t = zeros(height(data),1);
data(data.datequarter_num <EVENT_QUARTER,:).t = data(data.datequarter_num <EVENT_QUARTER,:).datequarter_num -  START_QUARTER;
data(data.datequarter_num >=EVENT_QUARTER,:).t = data(data.datequarter_num >=EVENT_QUARTER,:).datequarter_num -  EVENT_QUARTER;
data.IS_POST = data.post == 2;
nb_muz = length(unique(data.datequarter_num));

groupsummary(data,'datequarter_num','mean','assets')
groupsummary(data,'datequarter_num')
data2 = groupsummary(data,'datequarter_num');

data3 = join(data(:, {'datequarter_num','post','rssdid'}), data2);

data2 = sortrows(data,6);
% groupsummary(data,'datequarter_num','mean','assets')

q_bar = 50;
qmin   = 0.001;



% loglik_f_power_with_trend(data,EVENT_QUARTER, q_bar, qmin, betaq, sig, k, mu0, mu1)
% loglik  = @(b) loglik_f_power_with_trend(data, EVENT_QUARTER, q_bar, qmin, b(1), b(2), b(3), b(4), b(5));
% loglik  = @(b) loglik_f_power_with_varying_betas(data, EVENT_QUARTER, q_bar, qmin, b(1:nb_muz),b(nb_muz+1), b(nb_muz+2), b(nb_muz+3));
% ( apost,apre, q_bar, qmin, betaq, sig, k)
betaq = 1.6;
c     = 0.6;
q_up  = 52;
mu0   = 0;
mu1   = 0;


A       =  [];
b       = [];
Aeq     = [];
beq     = [];
nonlcon = [];
options = optimset('Diagnostics','on','Display','iter', ...
                    'TolX',1e-30,'TolFun',1e-30,'MaxIter',1000,'Algorithm','interior-point','MaxFunEvals',1000);

% with k as the parameter
loglik  = @(b) loglik_f_power(apost, apre, q_bar, qmin, b(1), sigparam, b(2));
x0      = [betaq;0.01];
xstart  = x0;
lb      = [1e-9;0];
up      = [Inf;Inf ];
[theta,stderr,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% with q_up as the parameter
loglik_f  = @(apost,apre, q_bar, qmin, betaq, sig, q_up) [log(f_density_power(apre, q_bar, qmin, betaq, sig, q_bar));... 
    log(f_density_power(apost, q_bar, qmin, betaq, sig, q_up))];

loglik  = @(b) loglik_f(apost, apre, q_bar, qmin, b(1),sigparam, b(2));
x0      = [betaq;q_up];
xstart  = x0;
lb      = [1e-9;q_bar];
up      = [Inf;Inf ];
[theta2,stderr2,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);



% Export table
FID = fopen(fullfile('parameters_estimates_baseline_MLE_50B_powerlaw_preDF_placebo.tex'), 'w');
fprintf(FID, '\\begin{tabular*}{\\textwidth}{l@{\\extracolsep{\\fill}}lll@{}}		\\hline   &  & Estimated value & S.E. \\\\ \\hline ');
fprintf(FID, concatenate_text(3,{'$\\beta$  & Exponent of the power law distribution    & ',theta(1),'     & [',stderr(1),'] \\\\ '}));
%fprintf(FID, concatenate_text(3,{'$\\sigma$ &  Measurement error volatility (in \\%%)  & ',sigparam*100,' & [','-','] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\exp(\\overline{q})$  &  Assets of marginal bank (\\$ Billion) & ',theta2(2),'     & [',stderr2(2),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{' $\\tau$  & Cost of regulation (\\%% of profit)      & ',theta(2) * 100,'& [',stderr(2) * 100,'] \\\\ '}));
fprintf(FID, ' \\hline \\end{tabular*} ');




%% Table 11, Panel C) Placebo test Using Only Post-DF data ($10B)

sigparam = 4.258/100;
% Starting with correct values: -------------------------------------------
amin = 3;
amax = 40;
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apost = data(data.post == 2,:).assets;
apre  = [];

EVENT_QUARTER = 2010.5;
START_QUARTER = min(data.datequarter_num);
data.t = zeros(height(data),1);
data(data.datequarter_num <EVENT_QUARTER,:).t = data(data.datequarter_num <EVENT_QUARTER,:).datequarter_num -  START_QUARTER;
data(data.datequarter_num >=EVENT_QUARTER,:).t = data(data.datequarter_num >=EVENT_QUARTER,:).datequarter_num -  EVENT_QUARTER;
data.IS_POST = data.post == 2;
nb_muz = length(unique(data.datequarter_num));

groupsummary(data,'datequarter_num','mean','assets')
groupsummary(data,'datequarter_num')
data2 = groupsummary(data,'datequarter_num');

data3 = join(data(:, {'datequarter_num','post','rssdid'}), data2);

data2 = sortrows(data,6);
% groupsummary(data,'datequarter_num','mean','assets')

q_bar = 20;
qmin   = 0.001;




% loglik_f_power_with_trend(data,EVENT_QUARTER, q_bar, qmin, betaq, sig, k, mu0, mu1)
% loglik  = @(b) loglik_f_power_with_trend(data, EVENT_QUARTER, q_bar, qmin, b(1), b(2), b(3), b(4), b(5));
% loglik  = @(b) loglik_f_power_with_varying_betas(data, EVENT_QUARTER, q_bar, qmin, b(1:nb_muz),b(nb_muz+1), b(nb_muz+2), b(nb_muz+3));
% ( apost,apre, q_bar, qmin, betaq, sig, k)
betaq = 1.6;
c     = 0.6;
q_up  = 11;
mu0   = 0;
mu1   = 0;


A       =  [];
b       = [];
Aeq     = [];
beq     = [];
nonlcon = [];
options = optimset('Diagnostics','on','Display','iter', ...
                    'TolX',1e-30,'TolFun',1e-30,'MaxIter',1000,'Algorithm','interior-point','MaxFunEvals',1000);

% with k as the parameter
loglik  = @(b) loglik_f_power(apost, apre, q_bar, qmin, b(1),4.145/100, b(2));
x0      = [betaq;0.01];
xstart  = x0;
lb      = [1e-9;0];
up      = [Inf;Inf ];
[theta,stderr,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% with q_up as the parameter
loglik_f  = @(apost,apre, q_bar, qmin, betaq, sig, q_up) [log(f_density_power(apre, q_bar, qmin, betaq, sig, q_bar));... 
    log(f_density_power(apost, q_bar, qmin, betaq, sig, q_up))];

loglik  = @(b) loglik_f(apost, apre, q_bar, qmin, b(1),4.145/100, b(2));
x0      = [betaq;q_up];
xstart  = x0;
lb      = [1e-9;q_bar];
up      = [Inf;Inf ];
[theta2,stderr2,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);





% Export table
FID = fopen(fullfile('parameters_estimates_baseline_MLE_20B_powerlaw_postDF_placebo.tex'), 'w');
fprintf(FID, '\\begin{tabular*}{\\textwidth}{l@{\\extracolsep{\\fill}}lll@{}}		\\hline   &  & Estimated value & S.E. \\\\ \\hline ');
fprintf(FID, concatenate_text(3,{'$\\beta$  & Exponent of the power law distribution    & ',theta(1),'     & [',stderr(1),'] \\\\ '}));
%fprintf(FID, concatenate_text(3,{'$\\sigma$ &  Measurement error volatility (in \\%%)  & ',4.145,' & [','-','] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\exp(\\overline{q})$  &  Assets of marginal bank (\\$ Billion) & ',theta2(2),'     & [',stderr2(2),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{' $\\tau$  & Cost of regulation (\\%% of profit)      & ',theta(2) * 100,'& [',stderr(2) * 100,'] \\\\ '}));
fprintf(FID, ' \\hline \\end{tabular*} ');





%% Table 11, Panel D) Placebo test Using Only Post-DF data ($50B)

sigparam = 2.291/100;
% Starting with correct values: -------------------------------------------
amin = 30;
amax = Inf;
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apost = data(data.post == 2,:).assets;
apre  = [];

EVENT_QUARTER = 2010.5;
START_QUARTER = min(data.datequarter_num);
data.t = zeros(height(data),1);
data(data.datequarter_num <EVENT_QUARTER,:).t = data(data.datequarter_num <EVENT_QUARTER,:).datequarter_num -  START_QUARTER;
data(data.datequarter_num >=EVENT_QUARTER,:).t = data(data.datequarter_num >=EVENT_QUARTER,:).datequarter_num -  EVENT_QUARTER;
data.IS_POST = data.post == 2;
nb_muz = length(unique(data.datequarter_num));

groupsummary(data,'datequarter_num','mean','assets')
groupsummary(data,'datequarter_num')
data2 = groupsummary(data,'datequarter_num');

data3 = join(data(:, {'datequarter_num','post','rssdid'}), data2);

data2 = sortrows(data,6);
% groupsummary(data,'datequarter_num','mean','assets')

q_bar = 40;
qmin   = 0.001;


% loglik_f_power_with_trend(data,EVENT_QUARTER, q_bar, qmin, betaq, sig, k, mu0, mu1)
% loglik  = @(b) loglik_f_power_with_trend(data, EVENT_QUARTER, q_bar, qmin, b(1), b(2), b(3), b(4), b(5));
% loglik  = @(b) loglik_f_power_with_varying_betas(data, EVENT_QUARTER, q_bar, qmin, b(1:nb_muz),b(nb_muz+1), b(nb_muz+2), b(nb_muz+3));
% ( apost,apre, q_bar, qmin, betaq, sig, k)
betaq = 1.6;
c     = 0.6;
q_up  = 40;
mu0   = 0;
mu1   = 0;


A       =  [];
b       = [];
Aeq     = [];
beq     = [];
nonlcon = [];
options = optimset('Diagnostics','on','Display','iter', ...
                    'TolX',1e-30,'TolFun',1e-30,'MaxIter',1000,'Algorithm','interior-point','MaxFunEvals',1000);

% with k as the parameter

%loglik_f  = @(apost,apre, q_bar, qmin, betaq, sig, k) [log(f_density_power(apre, q_bar, qmin, betaq, sig, q_bar));... 
%    log(f_density_power(apost, q_bar, qmin, betaq, sig, exp(q_upperbar_proportional(k, q_bar ,[]))))];

%loglik  = @(b) loglik_f(apost, apre, q_bar, qmin, b(1),sigparam, b(2));
loglik  = @(b) loglik_f_power(apost, apre, q_bar, qmin, b(1),sigparam, b(2));


x0      = [betaq;0.0009];
xstart  = x0;
lb      = [1e-9;0];
up      = [Inf;Inf ];
[theta,stderr,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% with q_up as the parameter
loglik_f  = @(apost,apre, q_bar, qmin, betaq, sig, q_up) [log(f_density_power(apre, q_bar, qmin, betaq, sig, q_bar));... 
    log(f_density_power(apost, q_bar, qmin, betaq, sig, q_up))];

loglik  = @(b) loglik_f(apost, apre, q_bar, qmin, b(1),sigparam, b(2));
x0      = [betaq;exp(q_upperbar_proportional( 0.0009,40,[]))];
xstart  = x0;
lb      = [1e-9;q_bar];
up      = [Inf;Inf ];
[theta2,stderr2,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);



% Export table
FID = fopen(fullfile('parameters_estimates_baseline_MLE_40B_powerlaw_postDF_placebo.tex'), 'w');
fprintf(FID, '\\begin{tabular*}{\\textwidth}{l@{\\extracolsep{\\fill}}lll@{}}		\\hline   &  & Estimated value & S.E. \\\\ \\hline ');
fprintf(FID, concatenate_text(3,{'$\\beta$  & Exponent of the power law distribution    & ',theta(1),'     & [',stderr(1),'] \\\\ '}));
%fprintf(FID, concatenate_text(3,{'$\\sigma$ &  Measurement error volatility (in \\%%)  & ',2.224,' & [','-','] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\exp(\\overline{q})$  &  Assets of marginal bank (\\$ Billion) & ',theta2(2),'     & [',stderr2(2),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{' $\\tau$  & Cost of regulation (\\%% of profit)      & ',theta(2) * 100,'& [',stderr(2) * 100,'] \\\\ '}));
fprintf(FID, ' \\hline \\end{tabular*} ');

%% Table 10, Panel A): Estimate the parameters $10B assuming lognormal distribution
% Starting with correct values: -------------------------------------------
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
data.year   = floor(data.datequarter_num);
muz       = {};
sigz      = {};
muz.pre   = 1 + mean(log(data(data.post == 0,:).assets));
muz.post  = 1 + mean(log(data(data.post == 2,:).assets));
sigz.pre  = std(log(data(data.post == 0,:).assets));
sigz.post = std(log(data(data.post == 2,:).assets));
sigz  =  mean([sigz.pre,sigz.post]);

amin = 3;
amax = 40;
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0| data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
data.year   = floor(data.datequarter_num);
Abar = 10;
nb_muz = 1;
EVENT_QUARTER = 2010.5;
apre  = []; %data(data.post == 0,:).assets;%data(data.post == 0,:).assets; % [];%;
apost = data(data.post == 2,:).assets;

q_bar = log(Abar);
q_up  = log(11);


A       =  [];
b       = [];
Aeq     = [];
beq     = [];
nonlcon = [];
options = optimset('Diagnostics','on','Display','iter', ...
                    'TolX',1e-30,'TolFun',1e-30,'MaxIter',1000,'Algorithm','interior-point','MaxFunEvals',1000);


% with exp(q_up) as the parameter
loglik  = @(b) [ log(f_density_lognormal(log(apre),q_bar, b(1)  , b(2), b(3), q_bar));
                  log(f_density_lognormal(log(apost),q_bar, b(1), b(2), b(3), log(b(4))))];

x0      = [ones(nb_muz,1)*2; sigz; 0.05; exp(q_up)];
xstart  = [ones(nb_muz,1)*2; sigz; 0.05; exp(q_up)];
lb      = [-Inf*ones(nb_muz,1);1e-100;1e-100;exp(q_bar)];
up      = [Inf*ones(nb_muz,1) ;Inf   ;0.1   ;Inf];

[theta,stderr,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% mytable = table(theta .* [1;1;1;1;1e3],stderr .* [1;1;1;1;1e3],xstart .* [1;1;1;1;1e3], 'VariableNames',{'EstCoef','SE','StartCoef'}, ...
%                 'RowNames',{'muz_pre','muz_post','sigz','sig','k'});
% final_theta = theta;


% with q_up as the parameter
loglik  = @(b) [ log(f_density_lognormal(log(apre),q_bar, b(1)  , b(2), b(3), q_bar));
                  log(f_density_lognormal(log(apost),q_bar, b(1), b(2), b(3), q_upperbar_proportional( b(4), exp(q_bar) ,[])  ))];

x0      = [ones(nb_muz,1)*2; sigz; 0.05; 0.01];
xstart  = [ones(nb_muz,1)*2; sigz; 0.05; 0.01];
lb      = [-Inf*ones(nb_muz,1);1e-100;1e-100;0];
up      = [Inf*ones(nb_muz,1) ;Inf   ;0.1   ;Inf];

[theta2,stderr2,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);


% Export table
FID = fopen(fullfile('parameters_estimates_baseline_MLE_10B_lognormal_constant_mu.tex'), 'w');
fprintf(FID, '\\begin{tabular*}{\\textwidth}{l@{\\extracolsep{\\fill}}lll@{}}		\\hline   &  & Estimated value & S.E. \\\\ \\hline ');
fprintf(FID, concatenate_text(3,{'$\\mu_{q}$  & Average undistorted log(asset)   & ',mean(theta(1:nb_muz)),'     & [',mean(stderr(1:nb_muz)),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\sigma_q$ & Undistorted log(asset) std.  & ',theta(nb_muz+1),' & [',stderr(nb_muz+1),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\sigma$ &  Measurement error volatility (in \\%%)  & ',100*theta(nb_muz+2),' & [',stderr(nb_muz+2)*100,'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\exp(\\overline{q})$  &  Assets of marginal bank (\\$ Billion) & ',theta(nb_muz+3),'     & [',stderr(nb_muz+3),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{' $\\tau$  & Cost of regulation (\\%% of profit)      & ',theta2(nb_muz+3) * 100,'& [',stderr2(nb_muz+3) * 100,'] \\\\ '}));
fprintf(FID, ' \\hline \\end{tabular*} ');



%% Table 10, Panel B) Estimate the parameters $50 assuming lognormal distribution
% Starting with correct values: -------------------------------------------
amin = 40;
sig = 0.05;
amax = Inf;
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
data.year   = floor(data.datequarter_num);
Abar = 50;
nb_muz = 1;
EVENT_QUARTER = 2010.5;
apre  = []; %data(data.post == 0,:).assets;
apost = data(data.post == 2,:).assets;

q_bar = log(Abar);
q_up  = log(50);


A       =  [];
b       = [];
Aeq     = [];
beq     = [];
nonlcon = [];
options = optimset('Diagnostics','on','Display','iter', ...
                    'TolX',1e-30,'TolFun',1e-30,'MaxIter',1000,'Algorithm','interior-point','MaxFunEvals',1000);


% with exp(q_up) as the parameter
loglik  = @(b) [ log(f_density_lognormal(log(apre),q_bar, b(1)  , b(2), b(3), q_bar));
                  log(f_density_lognormal(log(apost),q_bar, b(1), b(2), b(3), log(b(4))))];

x0      = [ones(nb_muz,1)*2; sigz; sig; exp(q_bar)];
xstart  = [ones(nb_muz,1)*2; sigz; sig; exp(log(55))];
lb      = [-Inf*ones(nb_muz,1);1e-100;1e-100;exp(q_bar)];
up      = [Inf*ones(nb_muz,1) ;Inf   ;0.04   ;Inf];

[theta,stderr,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);

% mytable = table(theta .* [1;1;1;1;1e3],stderr .* [1;1;1;1;1e3],xstart .* [1;1;1;1;1e3], 'VariableNames',{'EstCoef','SE','StartCoef'}, ...
%                 'RowNames',{'muz_pre','muz_post','sigz','sig','k'});
% final_theta = theta;


% with q_up as the parameter
loglik  = @(b) [ log(f_density_lognormal(log(apre),q_bar, b(1)  , b(2), b(3), q_bar));
                  log(f_density_lognormal(log(apost),q_bar, b(1), b(2), b(3), q_upperbar_proportional( b(4), exp(q_bar) ,[])  ))];

x0      = [ones(nb_muz,1)*2; sigz; sig; 0.001];
xstart  = [ones(nb_muz,1)*2; sigz; sig; 0.001];
lb      = [-Inf*ones(nb_muz,1);1e-100;1e-100;0];
up      = [Inf*ones(nb_muz,1) ;Inf   ;0.04   ;Inf];

[theta2,stderr2,~,~] = Max_lik(loglik,xstart,'Hessian',A,b,Aeq,beq,lb,up,nonlcon,options); % [theta,~,~,~] = fmincon(@(b) -sum(loglik(b)),xstart,A,b,Aeq,beq,lb,up,nonlcon,options);


% Export table
FID = fopen(fullfile('parameters_estimates_baseline_MLE_50B_lognormal_constant_mu.tex'), 'w');
fprintf(FID, '\\begin{tabular*}{\\textwidth}{l@{\\extracolsep{\\fill}}lll@{}}		\\hline   &  & Estimated value & S.E. \\\\ \\hline ');
fprintf(FID, concatenate_text(3,{'$\\mu_{q}$  & Average undistorted log(asset)   & ',mean(theta(1:nb_muz)),'     & [',mean(stderr(1:nb_muz)),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\sigma_q$ & Undistorted log(asset) std.  & ',theta(nb_muz+1),' & [',stderr(nb_muz+1),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\sigma$ &  Measurement error volatility (in \\%%)  & ',100*theta(nb_muz+2),' & [',stderr(nb_muz+2)*100,'] \\\\ '}));
fprintf(FID, concatenate_text(3,{'$\\exp(\\overline{q})$  &  Assets of marginal bank (\\$ Billion) & ',theta(nb_muz+3),'     & [',stderr(nb_muz+3),'] \\\\ '}));
fprintf(FID, concatenate_text(3,{' $\\tau$  & Cost of regulation (\\%% of profit)      & ',theta2(nb_muz+3) * 100,'& [',stderr2(nb_muz+3) * 100,'] \\\\ '}));
fprintf(FID, ' \\hline \\end{tabular*} ');


